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ABSTRACT 



Aims. We study the evolution of circumstellar massive disks around the primary star of a binary system focusing on the computation of 
disk eccentricity. In particular, we concentrate on its dependence on the binary eccentricity. Self-gravity is included in our numerical 
simulations. Our standard model assumes a semimajor axis for the binary of 30 AU, the most probable value according to the present 
binary statistics. 

Methods. Two-dimensional hydrodynamical computations are performed with FARGO. Besides the dynamical standard method to 
determine disk eccentricities, we apply a morphological method which may allow a better comparison with observations. 
Results. Self-gravity leads to disks that, on average, have low eccentricity. Moreover, the orientation of the disk computed with 
the standard dynamical method always librates instead of circulating as in simulations without self-gravity. The disk eccentricity 
decreases with the binary eccentricity, a result found also in models without self-gravity. 

Conclusions. Disk self-gravity appears to be an important factor in determining the evolution of a massive disk in a binary system. 
High eccentric binaries are not necessarily a hostile environment for planetary accretion. 

Key words. Planetary systems: formation; Planetary systems: protoplanetary disks; Methods: numerical 



1. Introduction 

All stages of planet formation in binary star systems might 
be strongly influenced by the distinctive gravitational field . 
dMarzari & SchollL l2000t lOuintana et all l2002t Pf hebault et all 
I2004L I20061 12008: IPaard ekooper et all I2008I) The protoplane- 
tary disk around a primary star, for instance, is perturbed by the 
companion affecting its morphology and dynamical structure. 
Spiral waves develop at major resonances in the disk . Its shape 
is expected to be elliptic with varying eccentricity dKlev et all 
2008; IPaardekooper et all l2008t iKlev & Nelsonl [2008b- Both 
these features may alter the dust sedimentation process on the 
mid plane of the disk and the grain ac cumula t ion in to planetes- 
imals which might even be inhibited (Nelson (2000)). More ec- 
centric trajectories of the grains, dragged by the gas, increase the 
impact rate but, on the other hand, might lead to destructive col- 
lisions by increasing the relative velocity. A more turbulent disk 
might, however, favor accumulation of dust in bigger bodies due 
to gravitational instability on a faster timescale. 

Understanding the morphology of circumprimary gas disks 
is also crucial for the subsequent phase of planetesimal accumu- 
lation. Indeed, several recent studies have shown that this stage 
might be severely hampered, or even stopped by the combina- 
tion of secular perturbations from the companion star and fric- 
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tion with gas dThebault et alJ (l2006l I2008L l2009t) : IXie & Zhou! 
(2008)). This combined effect forces a differential phasing of 
planetesimal orbits according to sizes, which leads to high im- 
pact velocities for a large fraction of planetesimal-planetesimal 
collisions, which could i n some cases lead to an acc retion hostile 
dynamical environment. Paardekooper et alJ d2008l) have shown 
that these results, obtained assuming static axisymmetric gas 
disks, could be further amplified when letting the gas disk evolve 
and "feel" the binary perturbations. The most accretion-hostile 
cases were obtained when the gaseous disk reached high ec- 
centricities, with impact velocities almost twice as high as in 
the axisymmetric case, while low eccentricity disks gave results 
similar to that with a static disk There is in principle a possi- 
bily for an eccentric gas disc to lead to lower drag effects, that 
is if this eccentricity is close to th e forced secular eccentrici ty 
of the planetesimals (see Eq.17 of Paardekooper et al. (2008)). 
However, this theoretical behaviour is only possible for a very 
specific radial profile of the gas disk density and was moreover 
never observed in any test simulations, indicating that this case 
is probably only a marginal possibility (see detailed discussion 
in Sec. 6.1 of IPaardekooper et all (12008)). 

Also the alternative formation mechanism for giant planets 
by rapid disk gravitational instability see ms to b e affec ted by the 
companion perturbations. According to iBossI (|2006i) the pres- 
ence of the secondary star mi ght induce planet f ormation even if 



his results are questioned in lMaver et al. 



(2006). 
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Before performing time consuming simulations of planetes- 
imal accr etion in a binary system with a hybrid code l ike, fo r in- 
stance, in [ Kiev & Nelson (20071). IPaardekooper et al.l d2008l) and 
Marzari_elalj£2008), we first investigate the disk evolution de- 
pending on the parameters of the system. This paper is devoted 
to the influence of the binary eccentricity and self-gravity on the 
disk structure. 

In addition, the eccentricity of circumstellar disks in bi- 
naries may also become in the future an observable feature. 
Next generation of interferometers (like ALMA, Atacama Large 
Millimeter/submillimeter Array) aimed to obtain high angular 
resolution millimeter and sub-millimeter images might resolve 
the molecular gas and dust compone nts of disks by deriving val- 
ues for the disk size and ellipticity. iLim and Takakuwal {2006) 
have already performed with VLA (Very Large Array) imag- 
ing of the dusty disks in the LI 551 IRS 5 binary system with 
a resolution of about 5 AU while with ALMA it might be possi- 
ble to gain a factor from 2 to 5. The morphology of protostellar 
disks might be derived by a comparison with numerical simula- 
tions. The matching of observable features will allow to derive 
important constraints on physical disk parameters. An important 
parameter for comparing observations and simulations is the ec- 
centricity of a disk. Classical image treatment methods devel- 
oped for measuring the eccentricity and orientation of elliptical 
galaxies, once applied to circumstellar disks may yield different 
values compared to those usually obtained by numerical simula- 
tions. We will discuss this problem in Section 2. 

Investigating the effect of a stellar companion on the evolu- 
tion of a disk surrounding the primary, necessitates the explo- 
ration of a very large parameter space determined in particular 
by physical parameters for the disk and orbital parameters of 
the binary system. In a first step, we focus on the influence of 
the binary eccentricity using for all other parameters mean stan- 
dard values which are expected from observations. For the sim- 
ulations, we use t he lat est release of the hydrodynamical code 
FARGO dM asset! (|2000|)) which include s now in particular disk 
self-gravity (Baruteau & Massed (120081) ). Comparison with for- 
mer simulations allows the investigation of its influence on the 
formation of structures in the disk, on its eccentricity and orien- 
tation. 

The paper evolves in the following way. We first compare 
in section 2 the standard way to derive the disk eccentricity in 
a simulation with image treatment methods developed for ob- 
served disks. In Sec. 3 we recall the numerical FARGO model 
and give the parameters of our simulated systems. Sec. 4 is de- 
voted to the importance of self-gravity on the evolution of a disk. 
Sec. 5 is focused on the dependence of disk evolution on binary 
eccentricity. In section 6 we summarize and discuss our results. 

2. Dynamical and morphological disk eccentricities 

Under the perturbations of the binary companion, the disk 
around the primary star modifies its shape. The three major per- 
turbing forces acting on each gas particle, which are the gravi- 
tational attraction exerted by the companion star, the gas pres- 
sure and the gravitational attraction by the disk, act in syn- 
ergy and modify the trajectories of each individual gas par- 
ticle. In a disk surrounding an isolated single star, trajecto- 
ries are in a first approximation nearly circular when the self- 
gravitation of the disk is neglected. The gravitational poten- 
tial due to the disk and the companion star modify the tra- 
jectories which can be approximated by elliptical trajectories 
with varying eccentricities. Hydrodynamical simulations of a 
disk yield for each cell of the grid the velocity of the gas 



which can be used to compute an individual eccentricity at- 
tributed to each cell. The mass-weighted average of the eccen- 
tricities of all cells has been u s ed as an estimator for the eccen - 
tricity of the disk dKlev et all (l2008l)l Pierens & Nelson] d2007l) . 
IPaardekooper et al] (120081) ). We call it in the following the dy- 
namical disk eccentricity ej. It is defined by: 



(D 



where e, is the eccentricity of each cell i of the grid, assuming 
that the local position and velocity vectors uniquely define a 2- 
body Keplerian orbit, while m; is its mass computed from the 
local disk density. is the disk mass. In a similar way, a mass 
averaged perihelion of the disk md with respect to the orientation 
of the binary's apsidal line can be computed. We call these the 
dynamical elliptical elements of the disk. 

If the self-gravity of the disk is included in the model, the 
calculation of e,- and m-, must be improved. When a two-body 
Keplerian orbit is derived for each gas cell it must account for 
the gravitational attraction of all the other disk components. An 
approximate perturb ative model can be used in this case as in 
Marz ari et al. (2008). The disk is approximated as a sequence 
of uniform density rings starting from the inner radius and out- 
wards until the outer border of the disk is reached. The density of 
each ring is calculated by averaging the local disk density within 
the ring derived from the numerical simulation. The gravitational 
attraction of each ring on each indiv idual gas cell can be analyt- 
ically computed (Krogh et al] (119821) ) and, being a radial force, 
it can be added to the central force of the star. A Keplerian orbit 
for each gas cell is computed by adopting a slightly increased 
mass for the star which must be updated at each timestep since 
the mass distribution within the disk changes because of spiral 
waves and disk eccentricity. 

By comparing the computation of the orbital elements of 
each disk component with and without the above mentioned al- 
gorithm to account for self-gravity we find that for eccentricities 
larger than 0.1 the discrepancy between a rough 2-body esti- 
mate of <?,■ and mi derived with a constant stellar mass and that 
obtained with a varying stellar mass because of the disk self- 
gravity are not very different; so we neglect this effect in this 
paper. We can also safely neglect the effects of pressure forces 
in a first approximation since they do not introduce a large error 
in the computation of a reliable value of both e, and nr,-. In con- 
clusion, hereinafter we compute the eccentricity e, of each disk 
component by using a 2-body approximation with fixed stellar 
mass. 

An important question is how to relate e,/ to the disk eccen- 
tricity an observer would obtain applying an image treatment 
method. Eccentricities of elliptically shaped objects, in particu- 
lar of elliptical galaxies, on a digitized image are obtained basi- 
cally by two classical methods. One method uses all pixels which 
form an object. Computing first and second order momenta of 
the object yields eccentricity, semimajor axis and orientation of 
an ellipse approximating the object. The second method uses 
contour lines formed by pixels with about the same intensity. An 
ellipse is then fitted to the contour line. We applied the second 
method, which is implemented in MATLAB, where we force 
the focus of the ellipse to lay on the star. This method yields the 
morphological parameters e m , m m and a m of a disk. FigQ]illus- 
trates the morphological method. Using the gas densities at each 
cell in a simulation, we produce a synthetic CCD image. The 
pixel values correspond to gas densities. The morphological pa- 
rameters depend on the choice of the intensity values of pixels 
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Fig. 1. Example of outer disk contour identification by 
MATLAB. The value of p is set to 1CT 6 . An ellipse is then fit- 
ted to the contour to compute e m , m m and a m . The color scaling 
gives the logarithm of the gas density in normalized units. 

and, hence, of a gas density value po forming the contour line 
which corresponds to the border of the disk. 

The advantages of using morphological parameters are the 
following: 

- They are independent of the internal structure of the disk 

- They are close to observable quantities 

- The morphological semimajor axis a,„ can be used as a reli- 
able estimate of the disk size. 

It is noteworthy that the value e m may significantly differ from 
the average value of the e,'s computed for a = a m because at the 
outer border of the disk the distribution of the e,'s and ?zr,'s is 
sparse. In Figf2]we show the distribution of both <?,■ and ex, as a 
function of the cell semimajor axis a,-. Close to the star the val- 
ues of vii are very similar independently of the azimuthal angle. 
However, at the outer border both the distribution open up and 
the values e, and cr, depend on the azimuthal angle. In the bot- 
tom plot we show a case in which the values of ex, are distributed 
in the range [0, 2n]. For this reason e m may be significantly dif- 
ferent from < e,(a m ) >. 

The morphological method to compute disk eccentricities 
depends on the choice of po which is not an objective choice 
but very personal. This is, however, a very common problem for 
observers applying image treatment methods. In addition, in our 
simulations we find that the border of the disk is always well de- 
fined and it drops to small values very quickly due to the tidal 
perturbations of the star. This is shown in Figj3] where at about 
8 AU from the central star the disk density rapidly drops. This 
makes us confident that the method to compute the overall disk 
eccentricity e,„ and its orientation m m is not very sensitive to the 
minimum density value selected in the computations when this 
is below 10~ 6 (obtained after exploring a significant numer of 
test cases) which will be the value of po always adopted in this 
paper. A test computation of e m for different values of po ranging 
from 10" 5 to 1CT 7 gave a maximum variation of less than 7%. 

The drastic truncation of the disk (see Fig|3j makes us also 
confident that the shape derived from the morphological analysis 
might be retrieved in the future from observations even if the 
disk is optically thick. The shape of an optically thin disk can be 




Fig. 2. Distribution of <?, (top plot) and tjt,- (middle plot) over the 
disk as a function of a, for both the SG and non-SG cases. In the 
bottom plot the density distribution averaged over the sectors is 
shown. 



retrieved by observations in the millimeter range and the column 
density (our model superficial density) is indeed proportional to 
the flux. For optically thick disks this is not necessarely true, 
however the sharp truncation of the disk might be retrieved by 
inspecting the scattered light from the star. 



3. Numerical set up 

3. 1 . The hydrodynamical model 

To model the evol ution of the disk we use the hydrodynami- 
cal code FARGO dMassetl F2000h solving the Navier-Stokes and 
continuity equations on a 2-dimensional polar grid. The mesh 
center lies at the primary and the indirect terms are included in 
the potential calculation. Since we are not interested in the evo- 
lution of the binary star system but only in the structure of the 
disk under the binary perturbation after the tidal truncation of 
the disk, we fix the orbital parameters of the companion star. In 
the latest release of FARGO a Poisson equa tion solver based on 
the Fo urier method has been implemented (Baruteau & Masset 
(2008)) allowing to properly model the disk self-gravity (here- 
inafter SG). The impact of SG on the disk evolution is relevant 
when modelling highly perturbed massive disks. 

In all the simulations we use a grid extending from 0.5 to 
15 AU from the primary star: the number of rings is 256 while 



4 



F. Marzari, H. Scholl, P. Thebault and C. Baruteau: Disks in binaries 



that of the sectors is 512. The aspect ratio h - H/r is constant 
all over the disk and set to 0.05. The initial density profile de- 
clines as S = £or~'^ 2 and is smoothly Gaussian truncated at the 
inner and outer borders. At the inner border S is reduced by a 
factor 10 in betwen 0.55 and 0.5 AU, while at the outer bor- 
der of the disk, starting from 1 1 AU, the density i s decr eased to 
E/ = 10-% within 2 AU. Following IKlev et al l d2008h . to pre- 
vent numerical instabilities at the outer border of the disk caused 
by very low density values, we introduce a density floor equal to 
Y,f. Whenever the evolution leads to a density value lower than 
£/, E is reset to Ey. In a series of tests this has proven to be a 
good value preventing instabilities and not leading to an artifi- 
cial mass growth of the disk. 

We adopt a density at 1 AU equal to 2.5e-4 in normalized 
units giving an initial mass of the disk equal to 0.0AM s u rl . Our 
disk is more massive than that studied in lKlev & Nelson! {2008) 
whose total mass was about 0.0015M V „„. Since the results do not 
scale only with the size of the system but also with the mass of 
the disk, we expect that the results can be different. 

We used a non-reflecting boundary condition in all our cal- 
culations. It is aimed at removing as many reflections as possi- 
ble off the grid boundaries , while allowing mass t o flow through 
the inner and outer edges. IKlev & Nelson! (2008) used a differ- 
ent kind of boundary condition to model the disk in y Cephei 
forcing a viscous inwards-flow of a disk at equilibrium at the in- 
ner border. With the non-reflecting boundary condition we have 
a mass inwards-flow which is comparable to that computed by 
IKlev & Nelson! d2008l) in terms of mass fraction of the disk and 
we do not observe any strong outflow of disk material during 
the periastron phase through the inner border. A small elliptic 
central hole in the disk develops with time (see FigJTJ, a con- 
sequence of the eccentric orbits of the gas in the proximity of 
the star. Since our grid extends down to 0.5 AU, all gas parti- 
cles with pericenter q lower than 0.5 AU are expected to exit 
the inner border of the disk. Due to the collimated values of the 
perihelia m the outflow of the particles with q < 0.5 AU creates 
a region of low density with an elliptical shape centered on the 
star (see Appendix A). If the zrr/'s were random, then the inner 
hole would have been circular with a radius larger than 0.5. In 
support of this interpretation we find that the size of the elliptic 
hole in the simulations is well reproduced by modelling its outer 
shape with a Keplerian orbit with the pericenter at 0.5 AU and 
with eccentricity equal to the average eccentricity e,- of the local 
gas cells. In conclusion, with our boundary condition, the code 
produces reasonable results taking into account the truncation of 
the disk at the inner edge. Of course, setting the inner rim of the 
disk to 0.5 AU is an overestimate since most observative mod- 
els assume that disk truncation occurs near corotation with th e 
inner region cleared by the magnetosphere (IShu et al.l {1994)). 
In this context, the inner elliptic hole in our simulations can be 
considered a numerical artifact since fluid elements having their 
pericenter at the inner edge, which is larger than the physical 
one, are not allowed to return into the grid. To test the influence 
of this issue on the disk eccentricity we performed a test simu- 
lation of our standard case with the inner border at 0. 1 AU. This 
should be a realistic simulation even if it appears impossible to 
get it far in time due to the short timestep imposed by the CFL 
condition. After 4 binary orbits the disk eccentricities in the two 
simulations differ by less than 5 % and the two density distri- 
butions, shown in Fig|3]appear similar. This makes us confident 
that, in spite of the larger inner hole produced in the simulations 
with the inner border at 0.5 AU, the models are still reliable in 
terms of disk eccentricity computation. 
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Fig. 3. Density distribution for the simulation with the inned ra- 
dius Ro = 0.1AU after 4 binary orbits compared to that with 
Ro = 0.5AU at the same time. The distributions are similar. We 
include also the initial density distribution for Rq = 0.5 AU. 



3.2. Our standard binary system 

We select our representative binary system by inspecting 
the orbital element distribu tions of binary systems given in 
Duq uennov & Mayor! (Il99ll) . From their histograms we derive 
the most probable orbital elements and mass ratio for most wide 
binaries in our neighborhood. The semimajor axis is chosen to 
be at - 30AU and the eccentricity et, = 0.4. The mass of the 
primary star is set to M p = 1M Q while that of the secondary 
is M s = 0.4M o . The orbital period of the companion star is 
about 134 yrs. This should represent the most frequent config- 
uration where planet formation may occur in binary systems. At 
the same time it is a very good example of a highly perturbing 
configuration due to the large eccentricity and mass of the binary. 
We adopt in most of our simulations a value of kinematic viscos- 
ity of 10~ 5 (normalized units) which corres ponds, at about 5 AU 
of the disk, to an a value of about 2.5 x 10" 3 ((Shakura & Sunvaevl 
dHH). 

To study the dependence of the disk eccentricity on eb we 
have varied this parameter from 0.0 to 0.6 at a constant step of 
0. 1 . We have also run a case without viscosity (inviscid case) to 
compare with the viscous case. 



4. Effect of self-gravity 

In our simulations we include the effect of self-gravity. The typ- 
ical parameter adopted to measure the relevance of self-gravity 
is the Toomre parameter Q: 



Q 



hM p 



(2) 



where r is the radial distance from the primary star, whose mass 
is M p , and E is the surface density of the disk. This formulation 
holds for our models where the disk is isothermal and it has a 
constant aspect ratio h = H/r. In Fig|4]we show the value of Q 
averaged over several rings around 4 AU for three simulations 
with different binary eccentricity (e^ = 0, 0.4, 0.6, respectively). 
In all cases Q is below 15 indicating that indeed the disk gravity 
may play a role in shapin g the disk. In addition, as suggested by 
Paardeko oper et all (120081) . the parameter that might really mea- 
sure the relevance of self-gravity is h ■ Q which is well below 1 in 
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Fig. 4. Average value of the Toomre parameter Q computed at 
r = 4AU from the primary star for different models with — 
(red dots), et, = 0.4 (green dots) and ei, = 0.6 (blue dots). This 
last case shows large variations due to the strong disk perturba- 
tions when the companion star is at periapsis. 



our cases. The large fluctuations of Q observed in the case with 
eh = 0.6 are due to the larger perturbations of the companion star 
when it is at the pericenter of a highly eccentric orbit. 

Fig |6] illustrates the effect of SG on the evolution of the ec- 
centricity, orientation, size and mass loss of the disk over about 
90 binary revolutions. Both the dynamical and morphological 
eccentricities ej and e m are significantly smaller in the case with 
SG and almost constant in the considered time interval, with ej 
being systematically larger than e m . The case without SG shows 
large oscil lations in qualitat i vely similar to the behaviour de- 
scribed in Paardekoope r et alj d2008l) (excited case) even if the 
disk mass and the binary parame ters of our standa r d cas e are 
different from those adopted by iPaardekooper et alj d2008l) . e m 
is lower but still it has a wide wavy pattern. Spikes in the disk 
eccentricity appear whenever the companion star approaches the 
periastron and the disk is more excited. 

A tentative explanation of why the dynamical and morpho- 
logical eccentricities of a disk are lower with self-gravity is pos- 
sibly due to the tidal interactions between the inner and outer 
portions of the disk. The inner disk acts like a tidal bulge on the 
outer p art tending to circularize t he orbits of the outer fluid ele- 
ments dMurrav & Dermota d2000)). This triggers a chain effect 
by which each inner ring forces a lower eccentricity on the 
outer one until a steady eccentricity distribution is reached. 
This stationary distribution has a lower eccentricity also in 
the inner regions, as it can be seen from Figj2](upper panel). 
There is also a delay between the inner bulge and the outer one, 
as it can be seen in Figj2] (middle panel). The perihelia of the 
inner gas cells are different from those of the outer ones and the 
transition is around 3 AU. The inner and outer parts of the disk 
exchange angular momentum and energy and, as a consequence, 
there is a tendency towards circularization of the orbits. This 
may lead to an overall lower eccentricity of a self-gravitating 
disk under the perturbations of an outer star. This explanation 
related to the tidal bulge appears more reasonable by inspecting 
Fig(5] where the evolution of a lower mass disk is illustrated. In 
this case the initial mass of the disk is 1/10 respect to our stan- 
dard case (Q ~ 26 in the central disk at t — 10 4 yr). The disk 



eccentricity still approaches a low value of dynamical eccentric- 
ity very close to that of the standard case also shown in the figure 
as a reference case but on a longer timescale with large oscilla- 
tions. Self-gravity acts on a longer timescale and the damping is 
slow, possibly because of the lower mass in the tidal bulge, but it 
is still effective. We performed also three simulations with pro- 
gressively decreasing intial mass and we find that self-gravity is 
still important when the mass of the disk is 20 times smaller than 
the standard case (Q ~ 58) but it does not affect the disk evolu- 
tion when the initial mass drops below 50 times that of the stan- 
dard case (Q ~ 150). In this latest case the behaviour of the disk 
eccentricity is close to that of the standard case without self- 
gravity in terms of disk eccentricity ej. This means that indeed 
the parameter meas uring the relevance of se lf-gravity might be 
hQ as suggested in Paardekooper et al. (2008|). 
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Fig. 5. Variation of the disk eccentricity e c i for a disk with ini- 
tialy 1/10 of the mass of the standard case (0.004 M ) compared 
to the standard case (0.04 M ). 



An additional effect that might explain why self-gravity 
keeps the disk eccentricity to a lower value is related to the 
pericenter passage of the star. As discussed in Section 5 and 
shown in Figj7] during the close approach of the star to the 
disk at the pericenter strong waves are excited and the overall 
shape of the disk becom es highly ellipsoid al. This phenomenon 
was already observed in iKley et all Q2008). Self-gravity might 
damp perturbations on the disk as it appears in a different phys- 
ical scenario concerning the effects of moonlets on Saturn rings 
iLewis & Stewart! d2009l) . A faster damping of the pericenter per- 
turbations of the companion star possibly leads to a less excited 
disk. 

Also the evolution of the disk orientation is significantly dif- 
ferent when self-gravity is turned on. As illustrated in Fig|6]top 
right panel, md librates around n instead of circulating with a pe- 
riod of about 5000 yrs as in the case without SG. A similar thing 
is observed for the morphological elements except that vj m , in 
the case with SG, librates around n/3. We will see in a subse- 
quent section that the center of libration for m m depends on the 
binary eccentricity. 

The total mass of the disk has a significantly different evo- 
lution in the case with SG as compared to the case without SG 
(FiglSl bottom left panel). The mass loss in the case with SG 
is significantly reduced and it is almost linear compared to that 
without SG. This confirms that with SG the disk relaxes into 
an almost steady state where its overall dynamical properties 
are almost constant in time or adiabatically changing on a long 
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timescale. On the other hand, without SG the disk appears to 
be more excited with its mass and eccentricity distributions still 
evolving after 12000 yrs. 

The semimajor axis a,„ describing the disk shape is about 7.8 
AU, which is larger than the critica l semimajor axis of 6.5 AU 
for dynamical stability computed by Holman & Wiegertl dl999l) . 
However, this is the critical semimajor axis for long term sta- 
bility (1 Myr) of massless particles under the action of gravity 
alone. It is reasonable that a disk, which is also affected by pres- 
sure forces and viscosity may behave slightly differently. In ad- 
dition, our simulations last only 1 .2 x 10 4 yrs while the stability 
limit of lHolman & Wiegertl £1999) is defined over a timescale of 
1 Myr. 

In conclusion, SG leads to an almost stationary disk on a 
short timescale and its internal structure appears to be more com- 
pact and less eccentric than in the non-SG case. The orientation 
of the disk librates around a fixed value, a behaviour similar to 
the ev olution of massive bod ies under a perturbative frictional 
force dMarzari & Scholll d2000h ). 

5. Dependence of disk eccentricity on e b 

Possible sources of disk eccentricity in our configuration are 
the forced component of eccentricity excited by the companion 
star, mean motion resonances between the companion star and 
disk gas particles (which include Lindblad and corotation reso- 
nance s) and viscous overstabilitv dKatoldl978l) : lLa tter & Ogilvie 
d2006h . 

We cannot attribute the disk eccentricity to the 3:1 mean mo- 
tion resonance (inner eccentr ic Lindblad resonance with m=2), 
as described in Lubowl d 19911) . in all our cases. Our scenario in- 
cludes a massive (M s = O.4M ) and eccentric companion star 
which truncates the disk inside the 3:1 m e an mo tion resonance. 
As described in Artvmowicz and Lubow ( 1994), the disk trun- 
cation moves closer to the star for larger values of eb- A good 
measure of the expected truncation limit is given by the al- 
ready quoted 2-body stable zone derived bv lHolman & Wiegertl 
(1999). According to their numerical simulations, the limiting 
semimajor axis for massless bodies orbiting the primary star 
ranges from about 10.5 AU when e/, = 0.0 to 5.7 AU when 
eh = 0.4 and to 3 AU when et, = 0.6. These values reasonably 
reproduce the size of our truncated disks. Only in the case with 
eh — the outer edge of the disk may be marginally involved 
with this resonance. However, other resonances are within the 
disk and they are given in Table Q] with the values of the semi- 
major axes were they are located. These values are computed 
within the 3-body model without accounting for pressure forces 
that, however, do not dramatically change the location of the res- 
onances, at least for what it concerns the present reasoning. 

While we may be confident that for low orders i the loca- 
tion is close to the effective one, when the value of k = i — j 
becomes large, the pericenter frequency to may play an active 
role in moving the location of the resonance from that we have 
estimated. 

Concerning the resonance effect, one would expect a reduc- 
tion of the disk eccentricity for increasing eh since the number 
of resonances inside the disk is smaller. On the other hand, for 
larger values of eh resonances are stronger and, in spite of their 
large values, they might anyway cause some instability in the 
disk. 

We do not expect a signficant contribution from viscous 
overstability in exciting disturbances in the disk since in the stan- 
dard case we observe a similar behaviour for different values of 
viscosity. 
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Table 1. Location of the mean motion resonances within the cir- 
cumprimary disk up to order 10. The critical argument is of the 
type iA c - jA g - nm c - mzn g where A c is the true anomaly of the 
companion star and A g that of the gas while m g is the perihelion 
longitude of the gas. m c is the perihelion longitude of the star 
which is constant in our simulations. 



An additional effect particularly relevant when the binary ec- 
centricity eb is larger, is the strong gravitational perturbations 
excited by the companion star during the pericenter passage. In 
Fig. [7] we show the strongly altered shape of the disk when the 
companion star is passing at the pericenter of the binary orbit. 
This strong disturbance is slowly damped when the star departs 
from the pericenter. This effects seems somehow to suggest that 
a higher value of eh favors larger values for both ej and e m . 

On the other hand, a companion star in a very eccentric orbit 
spends significantly more time far away from the disk's vicin- 
ity as compared to the case when the companion moves on a 
circular orbit. The effects of the strong perturbations during the 
periastron passage are damped down on a longer timescale. As a 
consequence, it is not obvious to predict 'a priori' the effects of 
the pericenter passage of the compan ion star on the disk eccen- 
tricity for different values of e/,- 

In order to test the dependence of both ej and e m on eb we 
have performed 7 long term simulations with values of eb rang- 
ing from 0.0 to 0.6 with fixed values for ah, viscosity and initial 
disk mass. We included SG in all the models. In Fig|8]we show 
3 selected cases with eb = 0.0, e/, = 0.4 (our standard case) 
and eb = 0.6. Both the dynamical and morphological eccentric- 
ities ej and e m indicate that larger values of e/, lead, on average, 
to lower disk eccentricities. It implies that the stronger distur- 
bances induced by a more eccentric companion star during the 
fast perihelion passage are not sufficient to stir up on average the 
eccentricity of the disk. Eccentric binaries are less effective in 
exciting the disk eccentricity. 

The effect of eb on Wd is not significant and mj, librates 
around n for any value of eb- Analyzing the evolution of vj m , 
however, we notice that its behaviour strongly depends on e&. 
For eb = 0.0 and <?/, =0.1 the orientation of the disk described 
by m m circulates with a period of some hundred years, while for 
larger binary eccentricities it librates around a decreasing value. 
For eb = 0.4 it librates around n/2 while for e/, = 0.6 it librates 
approximately around n/3. Self-gravity forces the disk appear 
to behave similarly to a s mall body under the actio n of a dissi- 
pative force. As shown in Marzari & Scholl d2000l) . an increas- 
ing eccentricity of the perturber causes the periastron of a small 
planetesimal orbit, which is perturbed by gas drag, to be aligned 
with increasingly smaller values of m. 

The disk eccentricity measured by the morphological e m 
shows a behaviour similar to ej for different values of eb- 
However, e m is always lower on average than e^. The case with 



F. Marzari, H. Scholl, P. Thebault and C. Baruteau: Disks in binaries 



7 



eb = 0.0 shows large oscillations around 0. 1 with the same pe- 
riod of circulation of m m . For the case with eb = 0.6, points 
well above the average are observed. They are computed by the 
algorithm estimating e m when the disk is strongly perturbed in 
correspondence to the periastron passage of the companion star 
and they are not significant. 

The semimajor axis of the disk estimated by a m is decreas- 
ing as a function of et, as expected. We rec all here that the lim- 
iting s emimajor axis for stability derived by iHolman & W iegert 
d 19991) is smaller when compared to a m (see FigJS] bottom left 
panel) for any value of <?/,. The timescale we are considering here 
(1.8 x 10 4 yr, i.e. 130 binary revolutions) is possibly not long 
enough to fully destabilize the outer border of the disk. This is 
confirmed also by the time evolution of the disk mass. We are 
not expecting that the dis k border relaxes exactly at the limiting 
semimajor axis given by iHolman & Wiegertl (1 19991) because of 
the presence of pressure forces and viscosity which alter the dy- 
namics of individual gas molecules from a pure 3-body problem. 
However, it is a good reference value for the size of the disk. 

In Fig|8] bottom right panel, we illustrate the progressive 
mass loss of the disk. After the initial large shrinking of the disk 
because of the tidal truncation, the mass loss continues at differ- 
ent rates depending on eb- The responsible mechanisms are: 



- mass in-fall of disk material on the star caused by the vis- 
cous evolution. 

- mass stripping during the perihelion passage of the compan- 
ion star 

- progressive dynamical destabilization of the outer parts of 
the disk 
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Fig. 9. Evolution of the disk mass and eccentricity (e^) for the 
viscous and inviscid case. 



The most effective mechanism appears to be the viscous evo- 
lution. In Fig|9]top left panel we show the mass loss as a func- 
tion of time for our standard case and the corresponding invis- 
cid case. The inviscid model has a comparatively much slower 
mass loss which may be ascribed to the other two above men- 
tioned mechanisms. In addition, when the companion star passes 
through the pericenter, the eccentricity of the disk is temporarily 
excited and this may cause a further in-fall of material across 
the inner border because the gas particle perihelion gets closer 
to the star. 

The dependence of the disk eccentricity and size on the bi- 
nary eccentricity is summarized in FigJTfJfor all our simulations. 
Both ej and e m show a decreasing trend for increasing et, con- 
firming that the perturbations are stronger in the circular case 
where a tidal response is constantly forced. It is noteworthy that 
this behaviour is not typical of self-graviting disks only. We per- 
formed 3 simulations with eb = 0.0, 0.4 and 0.6, respectively, ex- 
cluding self-gravity from the model. The dynamical eccentricity 
ed was decreasing from an average value of 0.42 for eb = 0.0 to 
0.25 for eb = 0.4 and 0.2 for eb = 0.4. The size of the disk mea- 
sured by a m shows an almost linear reduction for larger values 
of eb- This is due to the increasing strength of higher order mean 
motion resonances of disk particles and the companion star. In 
FigHO] bottom panel, together with a m the location of all reso- 
nances up to order 10 are illustrated by dashed lines. The dif- 
ferent truncation radius related to eb also changes the response 
of the disk to eccentricity forcing of the companion. Although a 
larger <?/, leads to stronger forcing, this may be more than com- 
pensated for by the disk truncation radius being closer to the 
primary. 



6. Discussion and conclusions 

We have analysed the evolution of disks in binary star systems 
including self-gravity. We study their properties by using two 
different set of parameters. The first set has been widely used 
and computes the eccentricity and orientation of the disk as a 
mass weighted average over the disk of the individual cell or- 
bital parameters. This method outlines the dynamical properties 
of the disk and is termed by subscript d. The second set is more 
oriented to observers and it measures the morphological elliptic- 
ity and orientation of the disk. It is termed with subscript m and 
it is computed by numerically fitting the outer edge of the disk. 

An important ingredient in modelling the disk evolution is 
self-gravity. We show that it forces the disk to behave more like 
a solid body in terms of orientation of the disk. Instead of cir- 
culating, the disk's apsidal line librates around a fixed value in 
eccentric binary systems like a body moving under the gravita- 
tional pull of the two stars and a dissipative force. In addition, 
self-gravity does not allow the disk to become very eccentric. 
For a binary eccentricity of e/, = 0.4 both the dynamical and 
morphological eccentricity ej and e m of the disk are lower than 
0.1. This is an importan t result in terms of plane tesimal accre- 
tion since, as shown by Paar dekooper et all (120081) . lower values 
of the disk eccentricity give lower impact velocities than with an 
eccentric disk. This might lead to an environment which is less 
hostile to accretion, although, even in low ed cases, the relative 
impact velocity might be too high to allow protoplanet forma- 
tion. 

A somewhat unexpected result is the increase of disk ec- 
centricity with decreasing binary eccentricity, when keeping the 
same semimajor axis. The case in which the binaries are on 
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Fig. 10. Variation of the disk eccentricity (upper plot), measured 
by ed and e m , and of the disk size given by a m (lower plot) for 
different values of the binary eccentricity e/,. The dashed lines in 
the lower plot mark the location of the mean motion resonances 
with the companion star up to order 9 and degree 10. 



detail. In models with higher binary eccentricity, we obtain ellip- 
tically shaped empty regions which needs confirmation by more 
extended simulations which are out of the scope of this paper. 
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Appendix A: Origin of the internal elliptic hole of 
the disk 

When the gas particles in a disk are on elliptic orbits, imposing 
a circular inner edge to the disk may lead to the formation of a 
central elliptic zone of low density. This is an effect related to 
the Keplerian nature of the orbits. When the particles are at peri- 
center, they may pass through the inner edge of the grid and then 
they are lost. This is simulated in Fig. A.l left panel where we 
compute 10 s two-dimensional elliptical orbits (which can be as- 
sociated to gas particles) with semimajor axis ranging from 0.1 
to 3 AU. The eccentricity and pericenter longitude have similar 
values to those computed by the hydrodynamical code (e,, ?zr,) 
while the mean anomaly is selected randomly. During the or- 
bit generation, anytime the pericenter is lower than 0.5 AU (the 
inner limit of the grid we have used in the hydrodynamical sim- 
ulations) the particle is thrown away. The number density of the 
surviving particles is illustrated in Fig. A.l left panel and it is 
very similar to the outcome of FARGO shown in Fig. A.l right 
panel. Not only the shape of the elliptical low density region is 
similar, but there is also in both figures an over-density located 
beyond the apocenter of the elliptical shape. 



circular orbits appears to be the most perturbing configuration 
in terms of disk eccentricity and this outcome is typical also 
of non-self-gravitating disks. It is significantly larger than for 
eb = 0.4 and e/, = 0.6. This is an indication that the strong 
perturbations during the short timespan in which the compan- 
ion star is at the pericenter are damped down when the compan- 
ion is far away. The circular case is more effective in exciting 
a tidal response from the disk. This does not necessarily mean 
that highly eccentric binaries have a higher probability of form- 
ing planets. In fact, the direct gravitational perturbations of an 
eccentric companion on the planetesimals lead to larger relative 
velocities that may not be necessarily damped by a low eccentric 
disk. The outcome in this case strongly depends on the interplay 
between the drag force due to the disk and the forced component 
of th e orbital eccen t ricity of planetesimals due to the secondary 
star dThebault et alj d2006l) ). These results are a first step towards 
a more comprehensive study made with a hybrid code where the 
evolution of the gaseous component of the disk is computed with 
a hydrodynamical scheme including self-gravity while the plan- 
etesimal trajectories are computed with a N-Body integrator. 

The evolution of a circumstellar disk in eccentric binaries 
must be further explored by changing additional parameters like 
the mass and semimajor axis of the companion star. Also nu- 
merical issues should be investigated like the dependence of the 
results on the boundary conditions, even if the one we have 
adopted is widely used. The region near the inner edge of the 
disk, where the mass flows inwards, must be explored in more 



BH 



I 



Fig. A.l. Number density of a set of test particles in Keplerian 
orbits depleted of all those having pericenter lower than 0.5 AU 
(left plot) and outcome of FARGO (right plot). The color den- 
sity in the right plot gives the logarithm of the gas density in 
normalized units. 
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Fig. 6. Comparison between the evolution of a disk in a binary (a/, = 30 AU, eh = 0.4) with and without SG (Self Gravity). The plot 
on top to the left (a) shows the variation of e</ with time while that on the right (b) illustrates its orientation mj. The middle plots 
(c,d) give the corresponding behaviour of e m and m m , respectively. The bottom plot on the left (e) shows the mass of the disk vs. 
time while that on the right (f) illustrates the evolution of the disk size described by a m , the effective semimajor axis of the elliptical 
disk. In all plots the self-gravitating case is represented by green dots while that without self-gravity by red dots. 
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Fig. 7. Snapshot of the disk surface density for our standard case after 15000 yrs of evolution. The companion star is at the close 
approach and it excites strong spiral waves in the disk and it significantly alters its shape. 



12 



F. Marzari, H. Scholl, P. Thebault and C. Baruteau: Disks in binaries 



E 

CD 



(fl 

CO 




3000 6000 9000 12000 15000 18000 
Time (yr) 




3000 6000 9000 12000 15000 18000 
Time (yr) 




< 

E 
cn 




6000 9000 12000 15000 18000 
Time (yr) 




6000 9000 12000 15000 
Time (yr) 



18000 



Fig. 8. 



15 
14 
13 
12 
11 
10 
9 
8 
7 
6 
5 
4 

3000 6000 9000 12000 15000 18000 
Time (yr) 

Dynamical and morphological elements for different binary eccentricities 



e b =0.0 
e b =0.4 
e b =0.6 - 



I I I I I L 



3000 6000 9000 1 2000 1 5000 1 8000 
Time (yr) 



